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We study a family of nonholonomic mechanical systems. These systems 
consist of harmonic oscillators coupled through nonholonomic constraints. 
In particular, the family includes the so called contact oscillator, which has 
Q been used as a test problem for numerical methods for nonholonomic me- 

chanics. Furthermore, the systems under study constitute simple models for 
continuously variable transmission gearboxes. 

The main result is that each system in the family is integrable reversible 
with respect to the canonical reversibility map on the cotangent bundle. By 
using reversible Kolmogorov-Arnold-Moser theory, we then establish preser- 
vation of invariant tori for reversible perturbations. This result explains 
qq previous numerical observations, that some discretisations of the contact os- 

cillator have favourable structure preserving properties. 

CN 1 Introduction 

In this paper we study a family of continuous dynamical systems consisting of nonholo- 
nomically coupled oscillators. The main result is that the family is integrable and that 
• i-H integrability is stable under reversible perturbations. Although numerical discretisations 

are not discussed in the paper, our main motivation is to obtain a rigorous explanation of 
the numerical observations by McLachlan and Perlmutter [16] , that some discretisations 
for nonholonomic problems have structure preserving properties. 
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We continue the introduction with some background in Hamiltonian, nonholonomic, 
and reversible dynamics. In § 2 we give a presentation of the family of problems under 
study, and we state the main results. These results are then proved in § 3. 

The concept of Liouville-Arnold integrability is fundamental in the study of continuous 
Hamiltonian dynamical systems. Through Kolmogorov-Arnold-Moser (KAM) theory, 
Liouville-Arnold integrability implies a foliation of the phase space in invariant tori, not 
only for the system itself, but also for "nearby" Hamiltonian systems (see Arnold [2, 
Ch. 10] and references therein). 

By backward error analysis a numerical discretisation method for a continuous dy- 
namical system can, at least formally, be interpreted as the exact solution of a nearby 
system [10, Ch. IX]. Roughly speaking, the KAM result then states that symplectic 
methods applied to Liouville-Arnold integrable systems preserve perturbed invariant 
tori [20, 21]. This explains the excellent performance of symplectic methods for Hamil- 
tonian systems. 

In nonholonomic mechanics, there is no corresponding KAM theory, so the foliation in 
invariant tori are in general not stable with respect to perturbations (see the discussion 
in the end of [1, § 5.4.2]). Nevertheless, integrability of nonholonomic systems is a well 
established research topic. In the current literature there are three main techniques to 
obtain complete integrability of nonholonomic systems. One technique uses reduction 
by symmetry, and examples can be found in [5]. Next, the special case of systems of 
Chaplygin type is treated using a nonholonomic version of the Hamilton- Jacobi theory 
in [11, 7, 17]. Finally, one can appeal to the existence of an invariant measure, which 
is the technique used in [14] and [1, §5.4]. All these techniques depend on some special 
property of the nonholonomic system at hand, and our case is no exception. Indeed, the 
essential property we use is that the underlying differential equation is a non-autonomous 
linear system whose fundamental matrix belongs to SO(n). 

For discussions on numerical discretisations of nonholonomic problems we refer to [16, 
8, 9, 12, 13], to the monograph by Cortes Monforte [6, § 7] and to the references therein. 
Due to the lack of KAM theory for nonholonomic mechanics, it is unclear which structure 
to preserve when discretising nonholonomic systems in order for the modified differential 
equation to remain integrable. 

Let Q be a smooth configuration manifold of dimension n. The corresponding phase 
space is given by the cotangent bundle T*Q. Recall that T*Q comes with a canonical 
symplectic form, which, in local canonical coordinates {qi,Pi), is given by to = Yli^-Qi A 
dpi. In addition, T*Q comes with a canonical reversibility map (as does any vector 
bundle), namely the bundle map p: T*Q — > T*Q given in local coordinates by (qi,Pi) i— > 
~Pi)- A diffeomorphism ift of T*Q is said to be reversible if p o ij) = o p. A 
continuous dynamical system on T*Q with flow map is called reversible if the flow 
map is reversible. If X is the vector field that generates then the system is reversible 
if and only if Tp o X = —X o p, or equivalently, p*X = —X, where p*X is the push- 
forward of the vector field (see, e.g., [15, §4.3]). A reversible system is called integrable 
reversible if T*Q can be foliated into invariant tori, which are also invariant with respect 
to the reversibility map p. That is, if there exists a reversible change of coordinates into 
action-angle variables [10, Def. XI.1.1]. 
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There are various ways to define nonholonomic systems. For an overview, see the 
monograph by Bloch [4, Ch. 5]. In our setting, a nonholonomic system is defined by 
a triple (Q,D,H), where Q is a configuration manifold, D is a regular distribution on 
Q, i.e., a vector subbundle of TQ, and H is a Hamiltonian, i.e., a smooth function 
on T*Q. The constraints are given by q G CD. Equivalently, if D* C T*Q is the 
annihilator of D, then the constraints are {T>* , q) = 0. If the one- forms t\, . . . ,Tfc span the 
codistribution T>* , then the Hamiltonian version of the Lagrange-d'Alembert principle 
yields the governing equations (see [4, § 5.2]) 

dH 

dp; 

Pi = "% + Z^ A ^' (!) 

i 

where the quantities Tjj are defined by Tj = Yli T ijdqi, and Aj are Lagrange multipliers. 

We assume that the differential algebraic equation (1) defines an underlying ordinary 
differential equation (ODE) on the manifold 

M={(g,p)GT*Q;(D*,|^)=0}. (2) 

Notice that if H is regular, which we assume, then M is a fibre bundle over Q. M is a 
vector subbundle of T*Q if and only if H is quadratic in p (which is typically the case). 
It is closed under the canonical reversibility map p on T*Q if and only if the Hamiltonian 
is reversible, i.e., H op = H, and in this case the corresponding dynamical system on M 
is reversible, i.e., poi^r = ip^ t o p. We say that the nonholonomic system (1) is integrable 
reversible if its corresponding system on M is integrable reversible. 

2 Main results 

We consider a family of nonholonomic systems on the configuration space Q = M 3 , with 
a distribution of the form 

25 = { (qi,Q2, Q3, G ™ 3 ; f(q 3 ) qi + q 2 = } (3) 

where / is a smooth function on M, and a Hamiltonian function on T*1R 3 of the form 

Ho(q, P ) = \jz{— + k <d) + F(qs,P3), (4) 

i=l 
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Figure 1: Illustration of a continuously variable transmission gearbox for which the dy- 
namics is governed by the nonholonomic system (5). The belt between the 
two cones is translated along the shafts in accordance with the (93, /J3) subsys- 
tem, thus providing a varying transmission ratio. The belt is kept in a plane 
perpendicular to the shafts, so that the belt keeps a constant length. The 
variables q\ and q<i denote the angular deflections of the shafts. 



where m; > 0, h% > 0, e > are constants and F is a smooth, reversible and regular 
Hamiltonian function on T*M. From (1) we get the governing equations as 

Qi = — , Pi = -hqi + A/(g 3 ), 
mi 



<?2 = , P2 = -K2<72 + A, 

dF OF 

Q3 = P3 = , 
OP3 OQ3 



(5) 



= /(<73)<7i + Q2- 

From a mechanical point of view, this system consists of two linear oscillators, coupled 
through a single nonholonomic constraint which depends on the independent Hamilto- 
nian subsystem described by the Hamiltonian F(q^,p3). 

Remark 2.1. System (5) is a simple model for a continuously variable transmission 
(CVT) gearbox, see Figure 1 for an illustration. Indeed, the shafts are attached to 
spiral springs that are fixed to a chassis, and the system governing 53 is an arbitrary 
Hamiltonian system with Hamiltonian F, which explains the form of the Hamiltonian (4). 
Now the belt imposes a constraint between the rotation velocities q\ and q\ which, 
assuming that both cones have the same size, is determined by #3 as q^q\ = (1 — (73)92. 
If we assume that q% < 1 (which correspond to assuming that the gear ratio is finite), 
then we obtain a distribution of the form (3). As a result, the system (5) describes a 
CVT gearbox where the "driver" is continuously and periodically shifting according to 
the independent subsystem (93,773), and /((ft) describes the transmission ratio. 
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The conceptualisation of CVTs is attributed to Leonardo da Vinci in 1490. In 1879, 
Milton Reeves invented a CVT for saw milling, which he later used also for cars. Today, 
CVT gearboxes are used in a variety of vehicles, particularly small tractors, snowmobiles, 
and motorscooters. 

Theorem 2.2. If the (93,^3) subsystem described by the Hamiltonian F is integrable 
reversible, then the nonholonomic system (5) is integrable reversible, or, more precisely, 
there exists action-angle variables (a,b,c,9,ip) £ IR 3 x (M/Z) 2 for which the underlying 
ODE of (5) takes the form 

a = 0, 9 = u(a), 

6 = 0, £ = £(a), (6) 
c = 0, 

and for which the reversibility map is (a, b, c, 9, ip) 1— > (a, b, c, —9, —(f). 
The proof is given in § 3 below. 

Consider now a perturbed Hamiltonian H e = Hq + eG, where e > and G is a 
reversible smooth function on T*M 3 . We need a Lemma on perturbations of reversible 
nonholonomic systems, which we specialise to our case. 

Lemma 2.3. For small enough e, the underlying ODE of the nonholonomic system 
defined by (M 3 ,D,i7 £ ) is isomorphic to a reversible perturbation of the underlying ODE 
O /(R 3 ,£,F ). 

Proof. First, notice that H £ is reversible and for small enough e it is regular. Let 
M £ be the subbundle given by (2) with H = H £ , and let X £ be the corresponding 
reversible vector field on M e that generates the dynamics. For small enough e the map 
: M e — t- T> is a reversible bundle automorphism. Thus, X e induces a reversible vector 
field Xq e on Mo by 



V dp dp 

This vector field is an e perturbation of Xq, i.e., Xq^ = Xq + 0(e). □ 

An integrable system is called KAM stable if "almost every" invariant torus, in the 
sense made precise e.g. in [18, 19], is stable under either Hamiltonian or reversible per- 
turbations. By using Theorem 2.2 and Lemma 2.3 we can now use the reversible KAM 
theorem on the system (5) to show KAM stability under reversible perturbations. 

Corollary 2.4. If 

1. the (^3,^3) subsystem is integrable reversible, 

2. f and F are real analytic, and 

3. the frequency functions uj (a) and £(a) in (6) are linearly independent, i.e., there is 
no constant c£l such that w(-) = c£(-), 
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then the foliation into tori of the underlying ODE of (5) is KAM stable under reversible 
perturbations. 



Proof. If the (93,733) subsystem is integrable reversible it follows from Theorem 2.2 that 
the underlying ODE of system (5) is integrable reversible. If /, F are real analytic it 
follows from the proof of Theorem 2.2 (see §3 below) that u> and £ are real analytic. From 
Lemma 2.3 it follows that a reversible perturbation of (5) corresponds to a reversible 
perturbation of the underlying ODE of (5). The KAM theorem, as given in [19, Thm. 1], 
then states that the system (6) is KAM stable under reversible perturbations if and only 
if the image of the map a t-t (w(a), £(a)J does not lie on a straight line passing through 
origo, which is exactly the last assumption in Corollary 2.4. □ 

Remark 2.5. The system (5) with e = 0, ki = = 1, F{q^,p^) = p\/2 + (/J/2 and the 
coupling function / is defined as f{q$) '■= 93, was used by McLachlan and Perlmutter [16, 
§5.1] as a nonholonomic test problem (called the contact oscillator) for various numerical 
discretisations. The authors showed that this system has quasiperiodic orbits, and they 
gave numerical results indicating that quasiperiodicity is preserved by some reversible 
integrators. Moreover, the authors study reversible perturbations of the type described in 
Lemma 2.3, and show by numerical experiments that the same reversible methods seem 
to preserve the invariant tori foliation, when the perturbations are small. By backward 
error analysis, a reversible method can be interpreted (up to exponentially small terms) 
as the exact solution of a reversible perturbation of the original nonholonomic system. 
Therefore, Theorem 2.2 and Corollary 2.4 provide a rigorous explanation of the numerical 
observations by McLachlan and Perlmutter [16]. 



The proof involves the following steps: 

1. Derive the ordinary differential equation constituting the dynamical system on Mo. 

2. Derive action-angle variables. 

3. Show that the transformation to action- angle variables preserves reversibility. 
3.1 Ordinary differential equation 

Our aim is to reduce the constrained system (5) into an ordinary differential equation. 
From the governing equation (5) it follows that P2/m2 = —f(Q3)Pi/ m i, which implies 
that the energy can be written 



3 Proof of Theorem 2.2 



p 2 fcigf + k 2 ql . . 



(8) 



where p : 
«2(<?3) : = 




i+^/fe) 2 \-i/2 



. We also define 
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so 



By differentiating the identity p = ai (g 3 ) Pi we obtain 

dp = ai{q 3 )— + 02(93) — , (9) 
mi m 2 

P = 01(93) — + 02(93) — • (10) 
mi m2 

By replacing p\ and p 2 in the equation above, and by substituting for p in the equation 
for <ji and fa, we now obtain the following system. 

"1(93) 

91 = P 

mi 

02(93) 

92 = P 

m 2 

p = -fciai(g 3 )gi - fc 2 o 2 (93)92 (11) 

93 = ^—(93,P3) 

i>3 = -^—(93,P3) 
^93 

We have thus obtained an ordinary differential equation describing the system (5) in 
the variables {qi,q2,Q3,P,P3) £ M 3 x M 2 , which parametrize M (which we think of as 
a trivial fibre bundle over Q = M 3 ). Since the Hamiltonian is reversible it follows 
that the canonical reversibility map on T*M 3 restricts to the map p: (qi,q2,Q3,P,Ps) 
(9i) 92) 93> —p, —P3), and that system (11) is reversible with respect to this map. 

3.2 Integrability 

First, since we assumed that the Hamiltonian subsystem (93,^3) is integrable reversible, 
let a £ R denote the action variable and 9 G IR/Z the angle variable, and u(a) the cor- 
responding angular velocity, so that the canonical reversibility map (93, p 3 ) i-» (93, — p 3 ) 
becomes (a, 6) 1— > (a,— 9). The system (11), with obvious abuse notation can now be 
written 

oi(a, 9) 

91 = P 

mi 

a 2 (a, 9) 

92 = P 

m 2 (12) 

p = -kiai(a, 9)qi - /c 2 o 2 (a, 9)q 2 
a = 0, 
9 = u(a). 

Let so (3) denote the Lie algebra of skew symmetric 3x3 matrices. The structure of this 
system becomes evident if we introduce the vector variable u = (\/kimiqi, \/fc 2 m 2 g 2 ,p) 
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and the so (3) valued function 



B(a,i 











-a 2 {a, 



(13) 



Then, in the variables (u, a,fl)eR 3 xlx E/Z, system (12) takes the form 

« = B(a, 0)tt 

d = 0, = cj(a). 



(14) 



Our aim below is to show that equation (14) is integrable for each fixed value of a, 
which then proves that equation (11) is integrable. If co(a) = 0, then B(a, 9) is a constant 
rotation matrix, so the system has invariant tori for these values of a. If u>(a) ^ then 
we can rescale the time variable so that the angular velocity of 6 is 1. The system (14) 
can then be written 



u = A a (6)u 
= 1. 



(15) 



where 



A a (0):=B(a,0)Ma). (16) 
That this system is integrable is a consequence of the following more general result. 

Theorem 3.1. Consider a differential equation of the form 



u 



A(e)u 
i 



(17) 



where (u,9) G M n x M/Z are the dependent variables and A: M/Z — > so(n) is a smooth 
mapping. There is a smooth change of variables 



* : R n x R/Z 3 (u, 6) ^ (v(u, 6), 6) G R n x M/Z 

such that system (17) expressed in the new variables (v,9) takes the form 

v = Av 
= 1 

for a constant matrix A £ so(n), whose spectrum a (A) is such that 

a(A) C [— i7r, i7r] . 



(18) 



(19) 



(20) 
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Proof. One defines the flow map 3>(t) as the solution operator of the differential equation 
defined for all r G M by 

f =A H , W (21) 

with initial condition at r = — the initial time matters because the differential equation 
is not autonomous. This means that if it? is a solution of (21), then 

w{t) = §(t)w(0) Vt G K, (22) 

and vice versa. 

Since A(r) G 5o(n) for all r, the flow map after one period, i.e., $(1), is an element 
of SO(n). As a result, there exists a matrix A G so(n), with the restriction (20) on the 
spectrum, such that 

$(1) = exp(A). (23) 
We define the mapping f: l n x I 1" x 1 by r) = (v(tt, r), r) and 

«(u,r) := exp(Ar)$(r) _1 M. (24) 

Consider a solution (tt(i),r(t)) of the differential equation 

tt(t) = A(r(t))u(t) 

r(0 = 1 (25) 

Define to = — T (0). Clearly we then have r(i) = t — to- By defining w(r) = u(t + to) 
we obtain w'(t) = u{t + to) = A(0(t + to)) = A(r), so to is a solution of (21). As a 
result, w(t) = <J>(t)io(0), which, using u(t) = to(r(t)), gives u(t) = <£>(t — to)u(to) = 
< I ) (r(t))u(to). We thus obtain that along a solution (it(t), r(t)) , 

$(r(t))" 1 w(t) = u(t )- (26) 

As a result, we obtain 

v(t) = Aexp(Ar)M(t ) = At>(t), (27) 

which proves the result for the mapping ^. 

Recall that, one of the pillars of Floquet's theory is that for any integer n G Z and 
any r G R we have ([3, § 28]) 

$(n + r) = $(r)$(ra). (28) 

Now, v(u, t) is periodic in r, of period one, because, due to the Floquet property (28), 
and the definition (23) of A, 

exp(A(r + l))$(r + l)- 1 = exp(Ar) exp(A)$(l) _1 $(T) _1 = exp(Ar)$(r)" 1 . (29) 

As a result, the mapping \t induces a mapping \I/ : R n xE/Z-> M n xl/Z which has the 
desired properties. □ 
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3.3 Reversibility 

We equip the space W 1 xl/Z with a linear involution R, defined as 

R(u,9) = (pu,-9), (30) 

for a given linear orthogonal involution p. One verifies by a straightforward calculation 
that the differential equation (17) is reversible with respect to R if and only if 

P A(e)p = -A(-e). (3i) 

Theorem 3.2. Assume that (31) holds, and that i— > A{9) is periodic of period one. 
Assume further that the spectrum <t(A) of the matrix A defined in Theorem 3.1 fulfils 

<t(A) C (-ivr,ivr). (32) 

ITien the mapping ^> : M ra x M/Z — >■ ]R n x M/Z, defined in Theorem 3.1, preserves re- 
versibility, i.e., Ro^f = ty o R. 

Proof. Using (31), one shows that 

p$(-r) = $(r)p. (33) 

This implies in particular that p<S>(— 1) = $(l)p, so using the Floquet property (28), we 
obtain 

p^(l)- 1 = *(l)p. (34) 
Now, recalling the definition of A in (23), we notice that (34) implies that 

exp(— A) = pexp(A)p = exp(pAp) (35) 

and, using the assumption on the spectrum of A, and the fact that the eigenvalues of 
pAp are the same as those of A, we deduce that —A = pAp. We therefore obtain 

exp(— Ar)p = pexp(Ar). (36) 

By combining (34) and (36) we get 

v(pu,—r) = exp(— At)<1?(— t)^ 1 pu 

= exp(-AT)p$(r) _1 u (37) 
= pexp(Ar)$(r) _1 M 

so we finally obtain 

v(pu, — r) = pv. (38) 
This finishes the proof. □ 

Remark 3.3. The condition (32) is a non-resonance assumption. We already know that 
the matrix A is chosen to satisfy (20), but we cannot exclude the extreme case of one or 
more half rotations. 
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3.4 Wrapping up 

We finish the proof by putting all the pieces together. Notice that with the canonical 
mapping p: M 3 — > M 3 defined as (qi, q2,p) (qi,Q2, —p), the property (31) is fulfilled for 
the matrix A a defined in (16). We thus obtain that, in the special case of R 3 , the map 
defined in Theorem 3.1 transforms the original system to an integrable one, and the 
map ^ is reversible. In M 3 , all rotations have one axis of rotation, so for all the values 
of the invariant a, we obtain two more invariants, and one angle variable. In the end, 
we thus obtain three invariants and two angle variables. 
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